home *** CD-ROM | disk | FTP | other *** search
/ The Very Best of Atari Inside / The Very Best of Atari Inside 1.iso / mint / mntinc20 / math-688.h < prev    next >
C/C++ Source or Header  |  1992-03-26  |  9KB  |  478 lines
  1. /******************************************************************\
  2. *                                   *
  3. *  <math-68881.h>        last modified: 18 May 1989.       *
  4. *                                   *
  5. *  Copyright (C) 1989 by Matthew Self.                   *
  6. *  You may freely distribute verbatim copies of this software       *
  7. *  provided that this copyright notice is retained in all copies.  *
  8. *  You may distribute modifications to this software under the     *
  9. *  conditions above if you also clearly note such modifications    *
  10. *  with their author and date.                               *
  11. *                                   *
  12. *  Note:  errno is not set to EDOM when domain errors occur for    *
  13. *  most of these functions.  Rather, it is assumed that the       *
  14. *  68881's OPERR exception will be enabled and handled           *
  15. *  appropriately by the    operating system.  Similarly, overflow       *
  16. *  and underflow do not set errno to ERANGE.               *
  17. *                                   *
  18. *  Send bugs to Matthew Self (self@bayes.arc.nasa.gov).           *
  19. *                                   *
  20. \******************************************************************/
  21.  
  22. /* Modified by Richard Stallman, November 1990, to initialize HUGE_VAL
  23.    specially on a Sun.
  24.    December 1989, add parens around `&' in pow.  */
  25.  
  26. #include <errno.h>
  27.  
  28. #ifndef HUGE_VAL
  29. #ifdef __sun__
  30. /* The Sun assembler fails to handle the hex constant in the usual defn.  */
  31. #define HUGE_VAL                            \
  32. ({                                    \
  33.   static union { int i[2]; double d; } u = { {0x7ff00000, 0} };        \
  34.   u.d;                                    \
  35. })
  36. #else
  37. #define HUGE_VAL                            \
  38. ({                                    \
  39.   double huge_val;                            \
  40.                                     \
  41.   __asm ("fmove%.d #0x7ff0000000000000,%0"    /* Infinity */        \
  42.      : "=f" (huge_val)                        \
  43.      : /* no inputs */);                        \
  44.   huge_val;                                \
  45. })
  46. #endif
  47. #endif
  48.  
  49. __inline static const double sin (double x)
  50. {
  51.   double value;
  52.  
  53.   __asm ("fsin%.x %1,%0"
  54.      : "=f" (value)
  55.      : "f" (x));
  56.   return value;
  57. }
  58.  
  59. __inline static const double cos (double x)
  60. {
  61.   double value;
  62.  
  63.   __asm ("fcos%.x %1,%0"
  64.      : "=f" (value)
  65.      : "f" (x));
  66.   return value;
  67. }
  68.  
  69. __inline static const double tan (double x)
  70. {
  71.   double value;
  72.  
  73.   __asm ("ftan%.x %1,%0"
  74.      : "=f" (value)
  75.      : "f" (x));
  76.   return value;
  77. }
  78.  
  79. __inline static const double asin (double x)
  80. {
  81.   double value;
  82.  
  83.   __asm ("fasin%.x %1,%0"
  84.      : "=f" (value)
  85.      : "f" (x));
  86.   return value;
  87. }
  88.  
  89. __inline static const double acos (double x)
  90. {
  91.   double value;
  92.  
  93.   __asm ("facos%.x %1,%0"
  94.      : "=f" (value)
  95.      : "f" (x));
  96.   return value;
  97. }
  98.  
  99. __inline static const double atan (double x)
  100. {
  101.   double value;
  102.  
  103.   __asm ("fatan%.x %1,%0"
  104.      : "=f" (value)
  105.      : "f" (x));
  106.   return value;
  107. }
  108.  
  109. __inline static const double atan2 (double y, double x)
  110. {
  111.   double pi, pi_over_2;
  112.  
  113.   __asm ("fmovecr%.x %#0,%0"        /* extended precision pi */
  114.      : "=f" (pi)
  115.      : /* no inputs */ );
  116.   __asm ("fscale%.b %#-1,%0"        /* no loss of accuracy */
  117.      : "=f" (pi_over_2)
  118.      : "0" (pi));
  119.   if (x > 0)
  120.     {
  121.       if (y > 0)
  122.     {
  123.       if (x > y)
  124.         return atan (y / x);
  125.       else
  126.         return pi_over_2 - atan (x / y);
  127.     }
  128.       else
  129.     {
  130.       if (x > -y)
  131.         return atan (y / x);
  132.       else
  133.         return - pi_over_2 - atan (x / y);
  134.     }
  135.     }
  136.   else
  137.     {
  138.       if (y > 0)
  139.     {
  140.       if (-x > y)
  141.         return pi + atan (y / x);
  142.       else
  143.         return pi_over_2 - atan (x / y);
  144.     }
  145.       else
  146.     {
  147.       if (-x > -y)
  148.         return - pi + atan (y / x);
  149.       else if (y < 0)
  150.         return - pi_over_2 - atan (x / y);
  151.       else
  152.         {
  153.           double value;
  154.  
  155.           __asm ("fmove%.d %#0rnan,%0"     /* quiet NaN */
  156.              : "=f" (value)
  157.              : /* no inputs */);
  158.           return value;
  159.         }
  160.     }
  161.     }
  162. }
  163.  
  164. __inline static const double sinh (double x)
  165. {
  166.   double value;
  167.  
  168.   __asm ("fsinh%.x %1,%0"
  169.      : "=f" (value)
  170.      : "f" (x));
  171.   return value;
  172. }
  173.  
  174. __inline static const double cosh (double x)
  175. {
  176.   double value;
  177.  
  178.   __asm ("fcosh%.x %1,%0"
  179.      : "=f" (value)
  180.      : "f" (x));
  181.   return value;
  182. }
  183.  
  184. __inline static const double tanh (double x)
  185. {
  186.   double value;
  187.  
  188.   __asm ("ftanh%.x %1,%0"
  189.      : "=f" (value)
  190.      : "f" (x));
  191.   return value;
  192. }
  193.  
  194. __inline static const double atanh (double x)
  195. {
  196.   double value;
  197.  
  198.   __asm ("fatanh%.x %1,%0"
  199.      : "=f" (value)
  200.      : "f" (x));
  201.   return value;
  202. }
  203.  
  204. __inline static const double exp (double x)
  205. {
  206.   double value;
  207.  
  208.   __asm ("fetox%.x %1,%0"
  209.      : "=f" (value)
  210.      : "f" (x));
  211.   return value;
  212. }
  213.  
  214. __inline static const double expm1 (double x)
  215. {
  216.   double value;
  217.  
  218.   __asm ("fetoxm1%.x %1,%0"
  219.      : "=f" (value)
  220.      : "f" (x));
  221.   return value;
  222. }
  223.  
  224. __inline static const double log (double x)
  225. {
  226.   double value;
  227.  
  228.   __asm ("flogn%.x %1,%0"
  229.      : "=f" (value)
  230.      : "f" (x));
  231.   return value;
  232. }
  233.  
  234. __inline static const double log1p (double x)
  235. {
  236.   double value;
  237.  
  238.   __asm ("flognp1%.x %1,%0"
  239.      : "=f" (value)
  240.      : "f" (x));
  241.   return value;
  242. }
  243.  
  244. __inline static const double log10 (double x)
  245. {
  246.   double value;
  247.  
  248.   __asm ("flog10%.x %1,%0"
  249.      : "=f" (value)
  250.      : "f" (x));
  251.   return value;
  252. }
  253.  
  254. __inline static const double sqrt (double x)
  255. {
  256.   double value;
  257.  
  258.   __asm ("fsqrt%.x %1,%0"
  259.      : "=f" (value)
  260.      : "f" (x));
  261.   return value;
  262. }
  263.  
  264. __inline static const double hypot (const double x, const double y)
  265. {
  266.     return sqrt(x*x + y*y);
  267. }
  268.  
  269.  
  270. __inline static const double pow (const double x, const double y)
  271. {
  272.   if (x > 0)
  273.     return exp (y * log (x));
  274.   else if (x == 0)
  275.     {
  276.       if (y > 0)
  277.     return 0.0;
  278.       else
  279.     {
  280.       double value;
  281.  
  282.       __asm ("fmove%.d %#0rnan,%0"        /* quiet NaN */
  283.          : "=f" (value)
  284.          : /* no inputs */);
  285.       return value;
  286.     }
  287.     }
  288.   else
  289.     {
  290.       double temp;
  291.  
  292.       __asm ("fintrz%.x %1,%0"
  293.          : "=f" (temp)            /* integer-valued float */
  294.          : "f" (y));
  295.       if (y == temp)
  296.         {
  297.       int i = (int) y;
  298.       
  299.       if ((i & 1) == 0)            /* even */
  300.         return exp (y * log (-x));
  301.       else
  302.         return - exp (y * log (-x));
  303.         }
  304.       else
  305.         {
  306.       double value;
  307.  
  308.       __asm ("fmove%.d %#0rnan,%0"        /* quiet NaN */
  309.          : "=f" (value)
  310.          : /* no inputs */);
  311.       return value;
  312.         }
  313.     }
  314. }
  315.  
  316. __inline static const double fabs (double x)
  317. {
  318.   double value;
  319.  
  320.   __asm ("fabs%.x %1,%0"
  321.      : "=f" (value)
  322.      : "f" (x));
  323.   return value;
  324. }
  325.  
  326. __inline static const double ceil (double x)
  327. {
  328.   int rounding_mode, round_up;
  329.   double value;
  330.  
  331.   __asm volatile ("fmove%.l fpcr,%0"
  332.           : "=dm" (rounding_mode)
  333.           : /* no inputs */ );
  334.   round_up = rounding_mode | 0x30;
  335.   __asm volatile ("fmove%.l %0,fpcr"
  336.           : /* no outputs */
  337.           : "dmi" (round_up));
  338.   __asm volatile ("fint%.x %1,%0"
  339.           : "=f" (value)
  340.           : "f" (x));
  341.   __asm volatile ("fmove%.l %0,fpcr"
  342.           : /* no outputs */
  343.           : "dmi" (rounding_mode));
  344.   return value;
  345. }
  346.  
  347. __inline static const double floor (double x)
  348. {
  349.   int rounding_mode, round_down;
  350.   double value;
  351.  
  352.   __asm volatile ("fmove%.l fpcr,%0"
  353.           : "=dm" (rounding_mode)
  354.           : /* no inputs */ );
  355.   round_down = (rounding_mode & ~0x10)
  356.         | 0x20;
  357.   __asm volatile ("fmove%.l %0,fpcr"
  358.           : /* no outputs */
  359.           : "dmi" (round_down));
  360.   __asm volatile ("fint%.x %1,%0"
  361.           : "=f" (value)
  362.           : "f" (x));
  363.   __asm volatile ("fmove%.l %0,fpcr"
  364.           : /* no outputs */
  365.           : "dmi" (rounding_mode));
  366.   return value;
  367. }
  368.  
  369. __inline static const double rint (double x)
  370. {
  371.   int rounding_mode, round_nearest;
  372.   double value;
  373.  
  374.   __asm volatile ("fmove%.l fpcr,%0"
  375.           : "=dm" (rounding_mode)
  376.           : /* no inputs */ );
  377.   round_nearest = rounding_mode & ~0x30;
  378.   __asm volatile ("fmove%.l %0,fpcr"
  379.           : /* no outputs */
  380.           : "dmi" (round_nearest));
  381.   __asm volatile ("fint%.x %1,%0"
  382.           : "=f" (value)
  383.           : "f" (x));
  384.   __asm volatile ("fmove%.l %0,fpcr"
  385.           : /* no outputs */
  386.           : "dmi" (rounding_mode));
  387.   return value;
  388. }
  389.  
  390. __inline static const double fmod (double x, double y)
  391. {
  392.   double value;
  393.  
  394.   __asm ("fmod%.x %2,%0"
  395.      : "=f" (value)
  396.      : "0" (x),
  397.        "f" (y));
  398.   return value;
  399. }
  400.  
  401. __inline static const double drem (double x, double y)
  402. {
  403.   double value;
  404.  
  405.   __asm ("frem%.x %2,%0"
  406.      : "=f" (value)
  407.      : "0" (x),
  408.        "f" (y));
  409.   return value;
  410. }
  411.  
  412. __inline static const double scalb (double x, int n)
  413. {
  414.   double value;
  415.  
  416.   __asm ("fscale%.l %2,%0"
  417.      : "=f" (value)
  418.      : "0" (x),
  419.        "dmi" (n));
  420.   return value;
  421. }
  422.  
  423. __inline static double logb (double x)
  424. {
  425.   double exponent;
  426.  
  427.   __asm ("fgetexp%.x %1,%0"
  428.      : "=f" (exponent)
  429.      : "f" (x));
  430.   return exponent;
  431. }
  432.  
  433. __inline static const double ldexp (double x, int n)
  434. {
  435.   double value;
  436.  
  437.   __asm ("fscale%.l %2,%0"
  438.      : "=f" (value)
  439.      : "0" (x),
  440.        "dmi" (n));
  441.   return value;
  442. }
  443.  
  444. __inline static double frexp (double x, int *exp)
  445. {
  446.   double float_exponent;
  447.   int int_exponent;
  448.   double mantissa;
  449.  
  450.   __asm ("fgetexp%.x %1,%0"
  451.      : "=f" (float_exponent)     /* integer-valued float */
  452.      : "f" (x));
  453.   int_exponent = (int) float_exponent;
  454.   __asm ("fgetman%.x %1,%0"
  455.      : "=f" (mantissa)        /* 1.0 <= mantissa < 2.0 */
  456.      : "f" (x));
  457.   if (mantissa != 0)
  458.     {
  459.       __asm ("fscale%.b %#-1,%0"
  460.          : "=f" (mantissa)        /* mantissa /= 2.0 */
  461.          : "0" (mantissa));
  462.       int_exponent += 1;
  463.     }
  464.   *exp = int_exponent;
  465.   return mantissa;
  466. }
  467.  
  468. __inline static double modf (double x, double *ip)
  469. {
  470.   double temp;
  471.  
  472.   __asm ("fintrz%.x %1,%0"
  473.      : "=f" (temp)            /* integer-valued float */
  474.      : "f" (x));
  475.   *ip = temp;
  476.   return x - temp;
  477. }
  478.